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O , Abstract. We present a detailed study of quantum simulations of coupled spin 

' systems in surface-electrode ion-trap arrays, and illustrate our findings with a proposed 

^H implementation of the hexagonal Kitaev model [A. Kitaev, Annals of Physics 321,2 

S (2006)]. The effective (pseudo)spin interactions making up such quantum simulators 

qh are found to be proportional to the dipole-dipole interaction between the trapped ions, 

'"^ and are mediated by motion which can be driven by state-dependent forces. The precise 

<y-\ forms of the trapping potentials and the interactions are derived in the presence of a 

J> surface electrode and a cover electrode. These results are the starting point to derive 

an optimized surface-electrode geometry for trapping ions in the desired honeycomb 
lattice of Kitaev's model, where we design the dipole-dipole interactions in a way 

^^ that allows for coupling all three bond types of the model simultaneously, without 

I- . the need for time discretization. Finally we propose a simple wire structure that can 

C^ be incorporated in a microfabricated chip to generate localized state-dependent forces 

^"^ which drive the couplings prescribed by this particular model; such a wire structure 

. . should be adaptable to many other situations. 

> 

X 

1. Introduction 

Ion trap systems have proven to perform weU for implementing the basic elements of 
traditional quantum computing, where evolution is described in terms of discrete gate 
operations, which can be implemented step by step as intermediate states are irrelevant. 
This is in contrast to quantum simulations, where the goal is to simulate the continuous 
evolution of a given Hamiltonian. While the initial proposal for quantum computing 
with trapped ions relied on a number of sequential steps to mediate effective qubit 
interactions [1], other approaches [2-10] achieve interaction between the internal states 
of the ions via constant Hamiltonians and therefore allow the development of quantum 
simulators based on trapped ions [7, 8, 11-15]. In such simulators, interactions between 
trapped ions are dominated by the Coulomb potential. For this interaction to affect 
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internal states (i.e., the qubits or pseudo-spins representing the effective quantum system 
to be simulated), state-dependent forces must be applied to some or all of the trapped 
ions. State-dependent forces can be achieved through optical ac Stark shifts [5, 7, 15-19], 
static magnetic-field gradients in combination with homogeneous radio-frequency (rf) 
fields [4, 10, 14, 20], or with rf field gradients [21, 22]. While in most cases the Coulomb 
interaction is considered between ions in a self-assembled single chain or crystal, coupling 
of independently trapped ions has recently been demonstrated [23, 24]. 

For quantum simulations with ions in microtraps, we must take into account how 
the presence of the electrodes modifies the Coulomb interaction. While in many systems 
this effect is negligible (for example, in the surface-electrode setup of [23] the Coulomb 
coupling was found to be enhanced by only 1.8 %, in agreement with our more general 
results in section 2.3), the general theory developed here for a lattice of surface-electrode 
(SE) microtraps shows that significant modifications to free-space couplings are possible. 
Far from being an inconvenience, these modified interactions can be used to design 
quantum simulations with specific short-range effective pseudo-spin interactions, which 
we illustrate with the hexagonal Kitaev model as a concrete example. 

The remainder of this paper is organized as follows. In section 2 we present a 
Green's function approach to solving electrostatic problems as they occur for surface- 
electrode ion traps in the presence of a cover electrode. In section 3 we derive general 
expressions for spin-spin couplings in two-dimensional microtrap arrays, applicable, for 
example, to electric coupling to light fields or magnetic coupling to microwave near- 
field gradients. In section 4 we combine all these methods to show how the hexagonal 
Kitaev model [25] can be implemented with an array of trapped ions on an optimized 
surface-electrode chip, including a dedicated wire structure that could be integrated in 
the chip to simultaneously mediate the couplings along three distinct bonds by use of 
magnetic-field gradients. Finally, Appendix A gives a summary of the used coordinate 
systems. 

2. Electrostatics in the presence of conducting planes 

The electrostatic interaction between charged particles close to conducting surfaces can 
be strongly modified by the presence of the conductors [27] . In the idealized geometry of 
a perfectly conducting grounded electrode plane aX z = 0, the total electrostatic energy 
of a set of charges Qi located at positions Vi in the half-space Zi > (see figure 1) is [26] 
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The Coulomb interaction term in (1) is expressed in terms of the Dirichlet Green's 
function Goo, which can be found from the free-space Green's function G^^\r^r') = 
l/||r — r'W by the method of images (see figure 1), 

G^{r, r') = ^ - ^ (2) 

V^p^ + {z- z'Y Vp2 + (^ + z'f 
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Figure 1. Coulomb interaction between two charges Q and Q' (full red circles) in 
the presence of a grounded plane. The image charges (empty red circles) are located 
below the grounded plane and carry opposite charge. Interactions between the charges 
(full blue arrow) contribute to (1) in full, while interactions between charges and image 
charges (dashed blue arrows) contribute with a prefactor of 1/2 [26]. 



where p = ^{x — x'y + {y — y'Y is the horizontal distance between the charges. 

In the following we review the effects of a grounded cover plane, i.e., a second, 
parallel conducting plane covering the electrode plane at height z = H (see figure 2a). 
In the initial proposal [28] and demonstration [29] of surface-electrode rf traps, the 
conducting surface nearest to the trap electrodes was theoretically at infinity but in 
practice a part of the surrounding apparatus. It has been suggested that adding a cover 
plane in the form of a dc-biased mesh above the electrodes could improve trap depth 
[30]. In addition to possible benefits of providing bias field and shielding, the cover 
plane could have more practical advantages, namely shielding the trapping region from 
fields due to quasi-static charges on insulators in the vacuum chamber, and establishing 
a more well-defined boundary condition. Further, if the cover plane is modified to 
carry rf and dc electrodes of arbitrary shape in the same way as the electrode plane, 
the presented formulas can be used directly to calculate the combined electric fields 
generated in this "sandwich trap" geometry (however, if optical access to such a trap 
geometry is achieved with holes and/or fiber optics in the electrode planes [31], the 
present full-plane treatment must be adapted [32]). 

Below, we first modify the Green's function (2) to include the cover plane and 
illustrate that a cover plane at height H leads to exponential shielding on a lateral 
length scale of H (section 2.1), then consider its effects on the electric potential generated 
by surface electrodes (section 2.2) and on effective dipole-dipole interactions between 
vibrating trapped ions (section 2.3). 



2.1. the shielding effect of the cover plane 

When a grounded conducting cover plane at height z = H is added to the setup of 
figure 1, the Coulomb interaction (1) of charges located between these two planes is 
modified to 
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Figure 2. (a) Sketch of a surface-electrode trap with a grounded cover plane 
positioned at a height H above the electrode plane. The red ring electrodes are at 
rf potential, while all grey areas are grounded. For static interactions or interactions 
varying slowly compared to the rf period, only the time-averaged potential contributes, 
so for our purposes the situation is equivalent to two completely grounded planes. 
(b) The interaction energy (3) between two point charges at same height h over the 
electrode plane, as a function of the charge separation p, in the presence of a cover plane 
at height H = lOO/i. The red and blue parts of the solid curves are computed by (4a) 
and (46), respectively, while the dashed lines illustrate the approximate behavior given 

by (5). 



Both the scaled self-potential cniz) and the Dirichlet Green's function Gh corresponding 
to the cover plane geometry with infinite conducting electrode planes at z = and z = H 
can be found by summing over an infinite sequence of mirror planes; and in the absence 
of a cover plane {H -^ oo) they reduce to (1) and (2). The scaled self-potential is 
euiz) = -[27 + ip{z/H) + ^(1 - z/H)]/{AH) = ^ + 0{zyH^) in terms of Euler's 
constant 7 = 0.577216 . . . and the digamma function tp{a) = r'(a)/r(a). The Dirichlet 
Green's function is 



GH{r,r')= J2 GUr + 2fiHz,r') 



(4a) 
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where Kq is the modified Bessel function of the second kind. The second form (46) 
is obtained by solving the Laplace equation in cylindrical coordinates [27]. Both forms 
converge for all parameters (p, 2;, z'); but while (4a) converges faster when Hr — r'|| < if, 
(46) is more suitable if ||r — r'|| > if, in particular for p ^ H, as discussed below. 

The Coulomb interaction energy Gh{p, z, z')QQ' /{Aneo) between two charged 
particles in a SE trap depends on the horizontal separation p, as illustrated in figure 2b. 
To illustrate this interaction energy we take the particles to be at the same height h 
above the electrode plane, and the cover plane height H to be much larger than h. 
When p <^ H, we expect the cover plane to be irrelevant, so that the interaction 
is described by a single image charge: it falls of as p~^ while p < h (where the 
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electrode plane is irrelevant) and as p~^ thereafter, as described by (2). When p ^ H 
the cover plane becomes important and the asymptotically dominant form is the first 
term of the resummation (46), so that the presence of the cover plane leads to an 
exponential shielding at the length scale of the cover plane height, as illustrated in 
Fig. 2. Summarizing, 

' 1/p for p ^ /i 

2/iVp^ for /K p < if 



GHip,h,h) s 
2.2. the potential due to the surface electrodes 



(5) 
sin2 /^^ j e-'^^/^ forp>iJ. 



The contribution to the total potential from the structured electrodes in the z = plane 
can be computed as an integral over the electrode plane: 

$(r)=/" G^H\r,r')<^ir')dx'dy', (6) 



where we have introduced a "surface Green's function" Gjj (r, r') = -^-^Gnir, r') | , . 
In the absence of a cover plane, the surface Green's function was found to be [32] 

G^^\ry) = G'i\p,z) = ^, (7) 

with the geometric interpretation that the potential at r due to an electrode at potential 
$0 is $o/27r times the solid angle spanned by the electrode as seen from r [33]. 
Alternatively, the electric field at r is proportional to the magnetic field that would 
be observed if a current were flowing along the edge of the electrode [33, 34]. For 
electrode configurations that are translationally invariant in the x-direction, the system 
can be described by conformally mapping the upper-half yz-plane {z > 0) to a disc [35]. 
Analogous to (4a) and (46), we have for the general case including a cover plane, 

oo 

Gf{p,z)= J2 G^£\p,z + 2^^H) (8a) 

/!= — oo 

-i.i:^-(^)-o(7) (8.) 
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where ( is the Riemann zeta function, Pj are Legendre polynomials, and s = [|r — r'|| = 
y^p2 + z^. Forms (8a) and {8b) converge for all {p,z)] but while (8a) converges faster 
when s < iJ, {8b) is more suitable ii s > H. Form (8c) is restricted to s < 2H and is 
most useful for s <^ H . Similar to (5) we find the approximate behaviors 



^\jh)-^\^-^: 



ioT p <^ z 



__,,,(_)e^-/'' forp»H. 



Quantum simulation of the hexagonal Kitaev model with trapped ions 6 

where ip'{a) = r"{a)/r{a) — ip'^{a) is the first derivative of the digamma function. We 
conclude that the influence of any surface electrode is exponentially damped at distances 
larger than H, which is advantageous for the experimental construction of quasi-infinite 
surface microtrap lattices in that it reduces the influence of the inevitable electrode 
boundary: at any point further than H away from the edge of the electrode and cover 
plane, the trap will look as if the electrode were infinitely large. 

Since the surface Green's function only depends on the x and y coordinates through 
r — r', (6) is a folding integral (convolution) [32] and can be rewritten as a product of 
the Fourier-transformed quantities, ^{k^, ky, z) = G\j {k^-, ky, z)^{kx, ky, 0), with 

1 r°^ 

l>(A;,, ky, z) = — Hx, y, z)e''^'^^+''yyMx dy (10) 

and a similar expression for the Fourier-transformed Green's function. The latter is 
cylindrically symmetric {k = ^Jk^. + kT), 

and allows a rather intuitive interpretation. All solutions of the Laplace equation with 
horizontal wavevector {/c^, ky\ are of the form e^'y^^^+^vv) [aj^e^^^ + a_e^'^^); the Green's 
function (11) gives the unique solution which has unit amplitude on the electrode plane 
[G]:^ (/c, 0) = 1] and zero amplitude on the cover plane \G\i {k,H^ = 0]. Therefore (11) 
gives the unique extension of a unit-amplitude potential plane wave from the z = 
plane into the z > half-space which satisfies the boundary condition of vanishing 
amplitude on the cover plane. The fact that the momentum-space representation of the 
surface Green's function (11) can be written without infinite sums greatly simplifies the 
description of infinite lattices of surface-electrode microtraps [36]. 

2.3. dipole-dipole interactions between trapped ions 

Trapped- ion quantum simulators couple internal degrees of freedom of the ions (typically 
hyperfine states or metastable D-states) through a state-dependent coupling to shared 
vibrational degrees of freedom [1, 6, 8, 14, 20, 21] (see section 3). A crucial ingredient 
of these couplings is the precise nature of the Coulomb interactions between the ions. 
Here we address the details of this latter point, since it will determine how to construct 
a quantum simulator of a desired system, as exemplified in section 4. 

We consider the regime of "stiff" ion trapping [6], where the Coulomb interaction 
is relatively small compared to the trapping potential, and we can interpret the normal- 
mode dynamics of the ion crystal as that of a set of local harmonic oscillators that are 
weakly coupled. The ion trapping potential defines a set of local eigenmodes for the i^^ 
ion corresponding to vibration in three orthogonal directions mf (with ||?Tif || = 1 for 
fi = 1,2, 3) around an equilibrium position Roi. In what follows we use these directions 
to parametrize the position of the i^^ ion as 

3 
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The total Coulomb energy of a set of N charges is given in (3), and the leading-order 
term that couples the motion of the ions is 
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Since we are mainly interested in near(est)-neighbor interactions, we evaluate this 
expression in terms of the infinite sum over image charge pairs (4a), rather than the 
resummed form (4&): 
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where the explicit dipole-dipole coupling is given by the expression without a cover 
plane, 

m ■ m' — 3(m. ■ n){m' ■ n) m ■ m' — 3{m ■ n){ffi' ■ n) 
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"')/||r — r'll, and the mirrored quantities 
f' = r' — 2{r' ■ z)z and m' = m' — 2{m' ■ z)z. The first term of (15) is the well-known 
dipole-dipole interaction, while the second term is the correction due to image charges 
in the electrode. 

In order to illustrate the behavior of the dipolar interaction (15) in close proximity 
of a conducting electrode plane, we again consider two ions located at equal height h 
above the electrode plane, spaced by a distance p along the x axis, and in the absence of 
a cover plane. If we assume that both ions vibrate along axes m = m' that are parallel 
to the lab-frame coordinate axes, then we find that the presence of the electrode plane 
can either increase or decrease the dipolar coupling strength: 
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Thus we see that by choosing the directions of vibration m in particular ways we can use 
the presence of the electrode plane to make the dipolar interactions fall off with the fifth 
power of distance instead of with the third power (which is the case in the absence of any 
conducting planes), as long as the ion oscillation frequency is low enough to avoid the 
effects of retardation and dissipation. The relevant length scale that determines whether 
or not the electrode plane has a strong influence on the dipole-dipole coupling is p ~ 2h, 
similar to figure 2; for even farther separations (p > H) we find exponentially damped 
dipole-dipole couplings due to the shielding effect of the cover plane (see section 2.1). 
These rapid dampings can be used to construct lattice simulation models with nearly 
local interactions, which is a desirable feature since many spin models from condensed- 
matter physics are formulated in terms of such local {e.g., nearest-neighbor) couplings. 
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3. Spin— spin interactions between trapped ions 

This section derives how state-dependent forces can induce pseudo-spin interactions 
between neighboring ions through the Coulomb potential. While this effect is well- 
known in principle [1], we show how these effective interactions are constructed in a 
lattice of ions without the need for time-slicing ( "Trotterization" [37]). Further we show 
that, to lowest order, the effective interaction strengths are proportional to the real- 
space Coulomb coupling strengths, an observation that greatly simplifies the design of 
lattice-based quantum simulators (see section 4). 

3.1. normal modes of vibration 

For small oscillation amplitudes rf the coupled harmonic motion of A^ ions in a lattice 
can be described by considering the local trapping potential curvatures (the second 
derivatives of the ion trapping pseudopotential with respect to position) around the 
equilibrium positions -Rqj and including the Coulomb couplings between ions in separate 
wells to second order [38]. If we assume that all excursions rj — Roi are already written 
in terms of the "bare" eigenmodes of the isolated local trapping potentials (with "bare" 
frequencies cuj^), as in (12), then the potential energy of the ions of mass M is 



V =-M 
2 



' N 3 N 3 



=1 /i=l i,j=^ Mi^=l 



(17) 



where 7,7 = and -f^j = ^^mf ■ ViVjGniRoi, Roj) ■ mT-, see (13). This quadratic 
potential energy can be diagonalized using coefficients Oi^m such that the real-space 
displacements can be written as 

3N 3N TV 3 

^i / ^ ^i^mQm WlXn ^ ^ ^iiim^jum ^ij^^u anCl ^ ^ / ^ ^i^m^i^rn' ^mm'y \^^) 

m=l m=l i=l IJ.=1 

and V = J2m=i k^^mlm i'^ terms of the lattice normal-mode amplitudes g™ and their 
frequencies Um- We quantize these normal modes through g™. ^-> ^m = <iom{,cim + ci\n) 
with gom = A//i/(2Ma;m) and the usual commutation relations [am,a|^,] = 5mm'- We 
will work in the "stiff" lattice limit, where we assume that the "bare" trap frequencies 
cjj^ = cD^ V? are equal for all ions| and the Coulomb couplings between ions do not 
significantly mix local modes with different values /i. This is the case when (i) the bare 
trap frequencies are sufficiently far apart: \'^ij\ ^ |w^ — w^| ^hJ,fJ' 7^ J^, and (ii) the 
vibrational bands are sufficiently narrow: {■Jifl <^ min,^ \cD^ — ul\ \/i,j,fj,. In this limit, 
the normal modes of the lattice separate into three disjoint sets indexed by ;U G {1, 2, 3}, 
each containing N normal modes with frequencies close to the corresponding u^. 

I This equality can be relaxed to the condition that for the modes that are used for inducing spin-spin 
couplings, the spreads of the bare frequencies are much smaller than the dominant Coulomb couplings. 
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3.2. state- dependent forces 

In addition to the Coulomb couplings, which are always on and define the coupled 
vibrational eigenmodes of the trapped ions, we can experimentally introduce fields that 
couple to internal states of the ions. Examples of such interactions are electric or 
magnetic dipole couplings, Raman couplings, or electric quadrupole couplings to laser 
or microwave fields. In the following general treatment we assume that there is a coupling 
between (a) classical external field(s) effectively oscillating with angular frequency ui, 
and two internal states of each ion, forming an effective two-level (spin-1/2 or pseudo- 
spin-1/2) system. Irrespective of the type of induced coupling, the coupling operator of 
the i^^ ion in its (pseudo)spin-l/2 subspace can be expressed as a linear combination 
of the identity operator (Tq and the Pauli matrices a) , i E {X,Y,Z} expressed in a 
quantization coordinate frame whose axes are given by the orthonormal vectors X, Y, 
Z (see Appendix A). The coupling Hamiltonian can thus be very generally expressed as 

N 

Ui^Y^ Y, [cf cos(c^it+0«) + (r,-i?o.)-sf cos{uJit+4'^)]crf\ (19) 
i=i te{o,x,y,z} 

where we have performed a first-order expansion in the ion positions assuming small 
oscillation amplitudes. Any type of spin-1/2 coupling that is used with trapped ions 
(including effective couplings to pseudo-spin degrees of freedom) can be brought into 
this form, where terms with non-vanishing prefactors c^ and sj, are referred to as 
"carrier" and "sideband" terms, respectively. The phases can absorb differences in the 
details of driving fields: while stationary fields in general have 0c = 0s , travelling 
waves {e.g., light fields) are characterized by 0c = 0s ± f • 

As an example, the coupling of physical spins to a magnetic field is found by 
expanding their magnetic dipole operators as ft^"^' = (?'^*^/iB((3"^ X -|- ayY + o^Z) in 
terms of the Bohr magneton /xb and the (yf-factors g^'^'^; for small ion excursions the 
coupling Hamiltonian to the magnetic field T-Lj = — ^^ p,"^' ■ B{ri) cos(wit -|- 0) can thus 
be expressed in the form of (19) with 
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» -0 



-(^(^VbX ■ B{R, 
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-9''V'BV|X-B(fi„i)| 



£ = 




£ = X, and similarly for l = 


--Y,Z 


i = 


(20) 

-y:z' 


i = X, and similarly for i = 



and 0c = 0s = Vi. We stress, however, that the above form of the magnetic dipole 
operator does not apply to pseudo-spins for their effective interactions with external 
fields; see section 4.4 for an example involving pseudo-spins. 
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3.3. effective spin-spin interactions 

Inserting the lattice normal-mode expansion (18) and (12) into (19), we can write the 
interaction Hamiltonian as 



n. 
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where we have dropped the approximation symbol and introduced Vt 
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Sf . It is common to transform into the interaction picture to assess 
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the dynamics induced by such an interaction Hamiltonian. In this transformation, the 
field-free Hamiltonian "Ho = J2m=i ^-^ 
dependence of the operators in (21): 



^^U Si=i ^z leads to a time 
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o"y H-)- cr^' sin(a;-f-|t) + ay cos(a;||t). 



(22) 



The terms involving a^ and (Xy can lead either to spin flips without affecting the motion 
("carrier" -transitions, mediated by c^ax and Cy uy and resonant at the frequency 
difference u^i between the pseudo-spin states) or to interactions that couple spins and 
motion ("sideband" -transitions or M0lmer-S0rensen interactions [3], mediated by s^ctx 
and Sy (Ty and resonant around W|| ± w^); the latter will dominate if they are not 
driven too strongly and \uji — Wf^ ± w^l ^ |wi — uj^i\ for one of the signs in ±. Here we 
concentrate instead on a drive with frequency \ui — u^\ <^ oo^ <^ w^l, close to one of the 
three bare eigenfrequencies of the uncoupled ion sites (we have chosen /i = 3 without 
restricting generality). In this case we can neglect the terms in s^ and Sy as they are 
far off-resonant, and all c), by careful design of the experiment (see section 4.4). The 
interaction Hamiltonian in the interaction picture thus reduces to a coherent drive 

N 3N 



i=l m=l 



after a second rotating- wave approximation, with the detunings Sm = Wm — ^i- 
Equation (23) can be exactly integrated via a Magnus expansion [39, 40] to yield the 
unitary evolution operator 



Uf\t) = exp 



X exp 



^ 3^ 1 iS t 
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3N 
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with (hy 
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. The first exponent describes a set of time-dependent coherent 
displacements to all normal modes that can entangle the motion with the internal 
pseudo-spin states. The second exponent constitutes a phase that depends on pairs 
of spins and can be interpreted as a spin-spin interaction. For a faithful simulation of 
interacting spins it is desirable that (A) the first term should be as close as possible to 



(21) 
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the identity operator, in order to avoid populating vibrational excitations, and (B) the 
second term should provide sizable phases for desired inter- ion couplings with i ^ j, as 
these represent the spin-spin interactions. It can be shown from the expression above 
or by the use of a canonical transformation [6] that (A) can be approximately met as 
long as \fiim£\ ^ \Sm\ for all {i,m,i). 

The above restrictions do not limit the time scale for simulations, as long as one 
assumes that sufficiently strong couplings can be induced by lasers or microwave field 
gradients. However, the energy scale of nearest-neighbor Coulomb interactions also 
plays an important role in determining simulation time scales, but this dependence is 
hidden in the normal-mode coefficients Oj^m of (18). To illustrate this point, we assume 
that fiimo = for all normal modes m and sites i [see (20) for an example]; but what we 
show below also holds for more general cases. Assuming negligible displacements [see 
point (A) above] the unitary evolution operator thus simplifies to 



Ui{t) = exp 



exp 



■ N 
N 



3N 



rW;v(i), 



CT^'O"^' COS 



O Z O z COS 
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Orri. 



fort >sup|5„^|. (25) 



In the "stiff" lattice limit (see page 8) we choose the drive frequency ui close to one of 
the bare frequencies, say 0)3, such that the detunings Sm will be much smaller for normal 
modes in this set than for the other normal modes; consequently, the sum over modes 
m in (25) can be restricted to an "active" set of A^ modes clustered around 003. If we 
further choose the drive frequency such that ^3 = 0)3 — ui is much larger than the spread 
of the normal mode frequencies in the active group, the series expansion 
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together with the relations in the "active" group of normal modes 
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(26) 



(27a) 



CXJh) 



m=\ 



simplifies the unitary evolution (25) to 



L^i(t) ~ exp 



1^03^ 



4^25. 



'rn\ ■ sfr 



X exp 
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where we have further approximated gom ~ ^03 = y^Ji^^Mcu^) . The first term in (28) 
is a global phase; it is the second term that mediates an effective spin-spin coupling on 
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the lattice of ions. It can be interpreted as the evolution under an effective spin-spin 
coupling Hamiltonian due to the driving of local mode /U = 3, 

N 

nfs = Y: J?-f^f (29) 

where the effective spin-spin coupling coefficients are 



?037ff COS (t)l^{rrv^^){rn] ■ s^p' 
8/icu3(5| 






We conclude that to lowest order the strength of the spin-spin coupling between two 
ions is determined by the geometric overlaps (m^ ■ s^ ) and (m.^ ■ s^ ) of the "active" 
local modes of vibration with the direction of the state-dependent force, as well as by 
the real-space Coulomb coupling strength 7^^ between the ions moving along these local 
modes [see (13) and (17)]. The relative phases (pl^ of the driving forces can be used to 
modulate the coupling strengths. These observations are used in section 4 to construct 
a quantum simulator on a lattice of trapped ions. Equations (29) and (30) faithfully 
describe the evolution of the system under the following conditions: 

• The vibrational band structure of the trapped ions must consist of clearly distinct 
bands which can be addressed individually (see page 8 for the conditions for "stiff" 
trapping). 

• The state-dependent force must be driven at a small detuning 6 from one of these 
bands: 6 must be large enough such that (26) is valid for this band, but small 
enough such that the contributions to (30) from other bands are negligible. 

• The amplitude of the state-dependent force must be small enough such that it 
does not significantly excite ion vibrations. The condition |f2jm^| ^ \6m\ given 
above implies that the state-dependent forces must be weaker than the force scale 
F ~ h\6\/qQ for the addressed band. 

The above arguments can be made in a very similar fashion for uxu^ and ay uy 
interactions by considering the interaction Hamiltonian (21) with uj ~ u^^ ± o)^ in the 
appropriate basis |±) = (|t) ±e^^||))/v^, where the considered spin-spin interaction is 



X 



diagonal. The only slight complication can arise from carrier terms proportional to c 
and Cy that are detuned by roughly the motional eigenfrequencies. For detunings from 
the sidebands on the order of the dipole-dipole interactions and correspondingly small 
drive strengths, however, these carrier terms can be safely neglected. 

To summarize this section, we have considered general effective spin-spin 
interactions in the limit of "stiff" ion trapping. We have shown that even in a lattice, the 
spin-spin coupling strength of any two ions depends on the dipolar Coulomb coupling 
between these two ions. To avoid appreciable entanglement between (pseudo)spins and 
ion motion, the detunings of driving fields need to be larger than the couplings they 
induce. This latter finding agrees with other work on simulation with trapped ions [8]. 
At the same time, the detunings cannot be much larger than the couplings between 
nearest neighbors, which determine the finer structure of the normal-mode spectrum 
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around the frequencies of the uncoupled ("bare") motion of an ion tightly bound in 
one of the trapping wells (see section 4 for a concrete example). These requirements 
impose stringent bounds on the time scales necessary to perform simulations. For 
example, in [23] two ions at a distance of 40 /im exhibited an exchange splitting of 
approximately 3 kHz, barely sufficient to demonstrate a few energy exchanges before ion 
heating profoundly altered the motion. Simulations that need to progress adiabatically 
with respect to this exchange period will therefore be experimentally challenging and 
may require reducing anomalous heating below what was measured in [23]. 

4. Kitaev model 

As an example of how to use the results of sections 2 and 3 in the design of a quantum 
simulator, we construct an implementation of the hexagonal Kitaev model [25] with 
microtrapped ions. In its ideal form, this exactly solvable two-dimensional spin model 
has a topologically ordered ground state with anyonic excitations, which makes it 
extraordinarily interesting for study in a quantum simulator with individual access to 
the constituent degrees of freedom. 

4.I. model and implementation 

The hexagonal Kitaev model [25] has the Hamiltonian 

X— links y— links Z— links 

defined on a honeycomb lattice of spin-1/2 particles, where the lab-frame bond vectors 
refer to Fig. 3: 

Ax = 40, 1,0} Ay = 4^3,-1, 0}/2 A^ = 4-y3,-l,0}/2.(32) 

In this way the Hamiltonian (31) associates each real-space bond direction (32) with 
a spin quantization direction; however, it is important to keep in mind that the bond 
directions and the associated spin quantization directions are not a priori related (see 
Appendix A). The ions are located on two sublattices Lo and L,, as shown in figure 3. 
Neighboring ions are a distance d apart. 

The form of (31) is exactly that of (29) summed over three concurrent driving force 
fields. As these driving fields will be at very different frequencies, they can be applied 
simultaneously in order to drive the full Hamiltonian (31). What is therefore needed in 
order to implement the Kitaev model is a set of "bare" vibrational directions of the ions 
such that the couplings ^^^ ^ and therefore the effective spin-spin couplings (30), match 
the particular geometry of the three terms in (31). 

We choose the ion trapping height to be half of the inter-ion distance, h = d/2, and 
the orthonormal principal axes of vibration for ions on the two sublattices as 



m 



Y 



{0,2,v^}/V6 m^ = {0,-2,V2}/V6 

{V3,~1,V2}/Vq ml = {-V3,1,V2}/Vq 

{-V3,-1,V2}/VQ m^ = {V3,1,V2}/VQ. (33) 



Quantum simulation of the hexagonal Kitaev model with trapped ions 



14 




Figure 3. Dipole-dipole interactions (15) with the central (green) site due to the 
vibrational directions of (33) for the case of /x = X, expressed as percentages of the 
dominant couphng [equal to 1.99 x Q"^ /{A-Keod?), see (34); values below 1% of this are 
not shown]. Figures 4 and 5 show the vibrational band structure induced by these 
couplings. The couplings for Y [Z) are found by rotating this figure by 120° (240°) 
clockwise. Red and blue wires are described in section 4.4. 



This particular choice of axes of vibration has the property that dipole-dipole couplings 
of the sort of 7^^-^ (z. e., coupling the mf vibration of the ion at -Rqj with the m^- vibration 
of the ion at Roj) are strongly dominated by the nearest-neighbor couplings required by 
the Kitaev model (31), shown in figure 3 for /i = X. These couplings can be calculated 
from (15) (in the absence of a cover plane); the resulting nearest-neighbor terms in the 
dipole-dipole coupling part of the Coulomb potential (17) are 



K, 



E 



Q' 



fi,u&{X,Y,Z} 



Aireod^ 



52 - 3v^ , 

2A 



■'/i,i^ 



48 ^ 



■'/j.,iyj 



E 






(34) 



In addition, there are dipole-dipole couplings to neighbors that are further away and 
that turn out to be larger than the off-diagonal terms (/x 7^ v) in (34). The vibrational 
normal-mode band structure due to all of these dipole-dipole couplings is shown in 
figure 4, with the effective density of states shown in figure 5. It consists of two bands, 
in which neighboring ions oscillate in-phase (upper band) and out-of-phase (lower band), 
and whose small frequency spread is indicative of the dominance of the nearest-neighbor 
coupling over all other couplings. 

Many dipole-dipole couplings of the sort of r'^r", with /i 7^ z/ are nonzero in 
this configuration; however they do not lead to effective spin-spin couplings if the 
underlying trap frequencies along the directions m^ and m.^ are strongly off-resonant 
(see sections 3.3 and 4.2). Thus neglecting any 11 ^ v couplings, the effective spin-spin 
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Figure 4. Vibrational band structure due to the dipole-dipole interactions of 
vibrations along one of the sets of axes in (33) (/x = X, see figure 3), corresponding to 
the density of states shown in figure 5. In the lower band (left) neighboring ions move 
out of phase; in the upper band (right) they move in phase. The first Brillouin zone is 
drawn in black. Frequencies (colors) are given in units of wq (see figure 5). 




Figure 5. Density of states of the vibrational bands (figure 4) due to the dipole- 
dipole interactions of vibrations along the sets of axes in (33), shown in figure 3. 
Left: the six bands consisting of three off- resonant doublets (/i = X, Y, Z) with center 
frequencies (bare trap eigen-frequencies) split by the golden ratio (see section 4.2); 
a) = (ljx^yojz)^ and ujqy/u} = 0.02 (much larger than in a realistic experiment). 
Right: zoom of one of the doublets. The two bands detailed in figure 4 are clearly 
visible, separated by ~ 4a;o, and show the extent to which the couplings of figure 3 
are dominated by the desired nearest-neighbor couplings. The scale of the bands is 
wom = Q'^/{8TTeoCJf,Md^). 
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Hamiltonian that is constructed from fi = X Coulomb interactions is approximately 
- Hx/Jx = E a« [af^-^ + 0.05 (a^'^^^ + a^'^^^) + . . 



+ \ Y. a«[0.06(a5+^--^-) + a5+^--^-))+...j55) 



ieLoUL, 

where the first sum contains couplings between the different lattices while the second 
sum contains couplings within the lattices; numerical prefactors for the perturbing terms 
are used for brevity, as in figure 3. Here we have assumed for simplicity that all ions 
are simultaneously pushed by the same state-dependent force with equal phase. The 
Hamiltonians "Hy and l-iz are found from (36) through rotations by ±120°, and the 
total effective spin Hamiltonian is 

ri = — ■Jx'l'X — Jyii-Y — Jzii-Zi (36) 

where Jx, Jy-, and Jz are effective coupling constants containing the diagonal coupling 
strength, the physical prefactors, as well as the mechanisms used for achieving these 
effective spin-spin couplings (see section 4.4). The topic of whether or not this slightly 
perturbed Hamiltonian (36) exhibits the same interesting topological phases as the ideal 
Hamiltonian (31), at zero or finite [41, 42] temperature, is beyond the scope of this 
article. We mention, however, that if the perturbative terms of (36) will be deemed too 
strong, they can be reduced further by driving the different wires with different relative 
phases or amplitudes [see (30)]. 

The presented configuration of trapping height and vibrational axes nearly 
maximizes the desired dipole-dipole couplings at the same time as it nearly mimimizes 
all undesired couplings. By numerical optimization we can identify a configuration 
that performs a few percent better than (33), but we have not been able to obtain an 
analytical description of this configuration. 

J^.2. surface- electrode trap design 

To have maximally incommensurate vibrational frequencies along the normal mode 
axes (33) we choose them in the golden ratio Ux '■ ooy '■ ooz = 4'~^ '■ ^ '■ 4> with 
= (1 + v5)/2. We use the algorithm of Ref. [36] to find an rf surface-electrode 
pattern that will generate an infinite honeycomb lattice of exactly such microtraps, 
with the following constraints: 

• The unit cell of the electrode pattern is defined by the vectors a = d{\/3, 0, 0} and 
6 = 4^3/2,3/2,0}. 

• The ion positions within the unit cell define the sublattices Lo and L,: Rqo = 
d{0, 0, 1/2} and Rq, = d{VS, 1, 1/2}. 

• The gradient of the rf electric potential generated by the surface electrodes must 
vanish at the ion positions in order to have minima of the rf pseudo-potential. 
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Figure 6. Left: optimized rf (blue) and do (white) electrodes for the constraints 
of section 4.2, with no spurious traps. Dimcnsionlcss trap curvatures are k = 0.080. 
Right: optimized electrodes for a honeycomb lattice with out-of-plane quadrupole 
confinement, trapping height h — d/2, and cover plane ai H = 50d. Dimensionless 
trap curvatures are k = 0.102. Coordinates as in figure 3. 



• The principal axes of the second derivative tensors of the rf electric potential at 
the ion positions are aligned with the directions given in (33), with eigenvalues 
proportional to {(f)^^, 1, —0} in the m,^,, m,^,, and m,f, directions, respectively. 

• A cover plane is located at a height H = 50d. 

The resulting electrode pattern is shown in the left panel of figure 6. It generates 
microtraps at the desired positions with dimensionless curvatures [36] k, = 0.080 and no 
spurious additional microtraps. This is to be compared with a simple out-of-plane 
quadrupole honeycomb lattice geometry (k = 0.102) as in Ref. [36] (see figure 6, 
right panel), which can potentially be deformed during the experiment via dc electrode 
potentials into satisfying the above constraints. Such dc electrodes might be necessary in 
any experimental implementation in order to null micro-motion of the ions [43] induced 
by manufacturing inaccuracies and stray charges. 



4-3. trap depth and trap loading 

The depth of the microtrap lattices generated by the electrodes shown in figure 6 are 
rather shallow. For the honeycomb lattice, which is more easily analyzed due to its p6m 
symmetry, figure 7 (solid line) shows the ponderomotive pseudo-potential along a vertical 
axis through any microtrap, in units of Epp = Q'^U^^/ {AMQ'^jd'^) . For ^Be"*" ions {Q = +e, 
M = 9 u) trapped with [/rf = 50 V and n,i = 2tt x 200 MHz in a lattice of rf = 30 /im 
{h = d/2 = 15 /im and H = 50d = 1.5 mm), we have Epp = 4.7eV = 5.5 x IO^'/cb K 
and thus a trap depth (ion-loss barrier) of 0.00185ii^pp = 8.7 meV = lOlfce K. This small 
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Figure 7. Total potential on a vertical axis through any microtrap formed by the 
electrode of figure 6 (right panel). The ponderomotive pseudo-potential is drawn with 
a sohd line, in units of _Epp — Q'^U^i/^AMn'^iCp). Two levels of dc biasing (either 
through the rf electrode or applying the same bias potential to the dc electrodes and 
the cover plane; vertically offset to leave the trapping minimum unchanged) are shown, 
in units of V^p = Epp/Q. For a bias voltage V > 0.125ypp the hexagonal lattice 
microtraps are the only attractors (trapping zones) for cooled ions. 



adiabatic trap depth is likely further reduced by the breakdown of the pseudo-potential 
approximation near the trap barrier. In order to reliably load these microtraps we can 
make use of the cover plane (see section 2.1): applying a positive dc bias potential to the 
dc electrodes and the same potential to the cover plane adds an electrostatic potential 
that pushes the ions towards the rf electrode and increases the trapping well depth (see 
figure 7). Since this dc potential is equivalent to applying a negative dc bias potential 
to the rf electrode, its dc electric field at the ion trap sites vanishes (by construction 
of the rf electrode shape) and it thus does not induce micro-motion [43] . We find that 
applying a small dc bias potential of at least 0.125Epp/Q ^ 0.6 V is sufficient to make 
the desired lattice of microtraps the only minima of the total potential (dashed line in 
figure 7). By applying a stronger bias voltage, the resulting total potential (dotted line 
in figure 7) is deep enough to trap ions produced by photoionization directly from a hot 
atomic beam. This bias will simultaneously cause the traps to be shallower in the xy 
plane. 



4.4- wires for magnetic interaction 

As described in section 3, effective spin-spin interactions between ions require internal- 
state-dependent forces to be applied to the ions. In the present model we propose 
to embed parallel wires below the electrode plane, which generate local magnetic field 
gradients at the positions of the ions [see (20)]. A relatively simple periodic grid of two 
different types of wires, indicated in red and blue in figure 3, suffices to implement the 
spin-spin interactions along all three bond types of the Kitaev model. As explained 
in section 3, one can induce pairwise u^ ax ^^^ o-y cry interactions by currents at 
frequencies that are near-resonant to u^i ± ux and w-f-^ ± uy (M0lmer-S0rensen type 
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interactions [2]); the ctz u* interactions can be driven with currents that are near- 
resonant to uz (phase-gate type interactions [5]). Because the three band- manifolds 
are well separated in frequency (figure 5), the dynamics of the three bond types can be 
driven simultaneously by currents at separate frequencies that are mutually off- resonant. 
The geometry of the wires is determined by the condition that we need to suppress 
the magnetic field at the position of all ions, in order to have negligible carrier 
interactions c^ in (20), while maintaining useful field gradients that couple to their 
target vibrational directions (33) to drive spin-spin interactions on all three bonds 
simultaneously [21, 22]. The magnetic field of two infinite sets of infinitely long wires 
parallel to the y axis as in figure 3, with distance d^ = d\/3/2 between wires of equal 
color, is 

y^o/biuc =^si^hg|-fsini^ ^^j^^, *sinhg| + fsin^- 



^V3 cos ^- cosh ^ dVS cos ^ + cosh ^ 



where hq is the magnetic constant. The field over the blue wires {i.e., at the ion 
positions) vanishes at height h„ if the ratio of currents is 

^ = _tanh2^. (38) 

J red dy/S 

With this current ratio the magnetic-field gradient at the ion positions is 

/ 1 




\/B^{x = nd^, z = h^) = — — - — smh ——=^ 

3a^ rfV3 



(39) 

VI o; 

for any n G Z. As this magnetic-field gradient decreases rapidly with increasing distance 
/i„ to the ions, one should place the wires as close as possible to the ions. On the other 
hand, they should not interfere with the trap electrodes. As a reasonable compromise 
for the following estimates, we can assume that the wires are below the electrodes such 
that h^ = d„ = h\/?>. The actual h^ in an experiment will probably be dictated by 
constraints in the micro-fabrication. 

We choose the quantization axis of the pseudo-spins of the ions to coincide with its 
associated bond direction, Z = A.z/d; however, any other choice of Z will be equally 
valid, and the experimenter's choice may depend on the available quantization fields. 



For our choice, sz = ^z x gji^ ^^n .2'"° siiih ^ ^^^ for all ions on both sublattices. 



blue 

and with 7^-^ = 4^^ j^f^s ^^^24 \^^^ diagonal term of (34)] the interaction strength (30) 
becomes 

^ 7o^/(^o^ ■ sz){m^ ■ sz) ^ 7r^(52 - 3^2) [qozOf^Bf^oQlj^^J? ^.^^^-4 27r/i„ 
^ 16Ma;|5| 432 ATreoMhwzS^zd'^ ^^"^ dy/s' 

where g is the effective g-iactoT such that the energy difference between the |t) and 
I4-) pseudo-spin states in a weak constant magnetic field along the quantization axis is 
AE^l = hu^i — gfisBz- -^biuc i^ ^^e current amplitude in the blue wires at frequency 
^z + ^z with \6z\ ^ (jJz- With d = 30 /xm, g = 1, M = 9u, Q = +e, and 
cuz = 27rx 5 MHz, this coupling strength is Jz = 7.6kHz x [I^f^J /Af[5z/i2nk}lz)]'^. 
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To avoid sizable entanglement of the (pseudo) spins with the ion motion, we need to 
fulfill 



J^ 



hSz [5^/(27rkHz)]3 ^ ' 

The hexagonal Kitaev model features interesting gapped phases with anyonic excitations 
of the ground state for example if |Ja;| = \Jy\ < iJzlf^, in which case the coupling 
constant of the resulting effective Hamiltonian is Jeff = J^Jy /{16\ J z\^) < 17x1/256 
{see [25] for discussions of these phases and their emergence from (31)}. 

While the geometric prefactor of (40) depends on the details of the model, its 
functional dependences are expected to remain the same for a broad class of wire-driven 
coupled pseudo-spin models, in particular also for Jx and Jy of the same system. The 
exact form of the interactions along the X and Y bonds depends on the transition dipole 
matrix elements n^ = (t|/i|4,) which can have components along all spatial directions. § 
The component along the quantization axis Z is relevant for vr transitions (where the 
Z component of the total angular momentum F of the ion does not change, Aitlf = 0) 
while perpendicular components can be used for a^ transitions (during which the Z 
component of F changes by Amp = ±1). The coupling strengths for vr transitions are 
found by scaling (40) to be 



Jx/Y — Jz 



%X/Y[.t^d 



qoZQfJ'B 
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rr(^z)l2 ' 
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with currents of amplitude I^^^^ ^^^ at frequencies u^i ± {ux/v + ^x/v) that are 
required to drive M0lmer-S0rensen-type interactions [2]. For ^^-transitions analogous 
relations hold involving the projections of fi^ along (X ± iY)/\/2. We conclude 
that all interactions are similar in magnitude and that the current amplitudes can 
be used to tune the effective coupling strengths Jx, Jy, and Jz of the Kitaev model 
simulator (36). To drive all bonds simultaneously, a total of five alternating currents 
at different frequencies are necessary; using the largest allowed value of 1 in (41), the 
maximum rms current each of the blue wires [and red wires, see (38)] has to sustain is 
\/(^biue(^)) ~ \/572 X 0.36 A X [5/(27rkHz)]3/2. ^^ ^^^^^^^ however, that (40) and (42) 
depend very strongly on the vertical distance h^, and even a small reduction in /i„ can 
substantially reduce the required currents. 

Expression (40) seems to suggest that for quantum simulators built with the 
principles described here, decreasing the physical size of the ion-trap lattice (as given 
by the length scale d) will strongly increase the simulation speed, which is given by 
the effective dynamics of the particular quantum simulator but ultimately proportional 
to the spin-spin coupling strengths. However, a careful analysis of assumptions and 
constraints, carried out below for interactions along Z but equally valid for all other 

§ In the case of a real spin, fi^ — ^gfiB{X — lY) (see the example on page 9); but what follows also 
applies to more general pseudo-spin cases where we can have fx^- Z ^ 0, for example if the pseudo-spin 
states are hyperfine states of an ion. 
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effective interactions, disproves this observation. Firstly, if we assume tliat avoiding ion 
motion is a significant experimental constraint, a constant ratio Jz/i^^z) in (41) implies 
that the effective spin-spin coupling strength scales as 

rr('^z)l2/3 

J^ ^ 2 /T . (43) 

for a given ion species and given electrode shapes. Secondly, assuming the currents to be 
limited by heat dissipation, an upper bound on I^^^J must scale as rf^/^. And lastly, lower 
bounds for the trap frequency can be found in two ways: (i) to meet our assumption of 
stiff trapping, i.e., cjqz <^ ujz, we require coz ^ (i~^/^ x \Q\/\/87reoM; and (ii) to use 
the expansion (26) we require Uqz *^ \^z\, which, combined with the scaling of (43) for 
6z and with the above current scaling, implies a scaling of d~^ for the lower bound of 
uz- Together, these bounds imply that the maximal achievable coupling strength (43) 
scales only as d^^^^ or even d'^, depending on which of the frequency bounds is more 
stringent. Since current experimental setups are far from reaching these lowers bounds 
on Uz, miniaturization is expected for now to increase the coupling strength faster than 
these estimates; but the optimal dimension, where the ratio of simulation speed and 
heating rate (anomalous heating, scaling as d~'^ [44, 45]) is maximized, remains an open 
question. 

5. Conclusions 

We have discussed modifications to Coulomb potentials and interactions of trapped ions 
due to the presence of trap electrodes and cover planes. For plane geometries we have 
treated these modifications rigorously, using the method of image charges. We have 
found considerable deviations of the long-range behavior from that in free space when 
the relevant distances are of the order of ion-to-surface distances or larger. Moreover, we 
have developed a general approach to treating the effective spin-spin interactions of ions 
trapped in a multi-trap array in the stiff-trapping limit, where dipole-dipole interactions 
between nearest neighbors produce only small corrections to the bare normal modes 
of a given trap well. We have shown that effective coupling strengths, and therefore 
simulation timescales, are determined by the nearest-neighbor dipole-dipole couplings. 
As an illustration of the versatility and power of this stiff-trap-array approach, we have 
discussed a quantum simulation of the hexagonal Kitaev model. We have also addressed 
several practical challenges, including how the trap depth of the array may be improved 
so ions created from a thermal source with large kinetic energies can be trapped. 
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Appendix A. Summary of used coordinate systems 

In order to help distinguishing the various coordinate systems used in the text, we 
summarize them here. 

• The laboratory frame is spanned by the unit vectors 

i = {1,0,0}, j^ = {0,1,0}, £ = {0,0,1}. (A.l) 

Its orientation is shown in figures 1, 2, and 3. In section 2 lab-frame vectors are 
written as r = xx + yy + zz with p = y^x^ + y^ and r = ^yx'^ + y"^ + z^. 

• The pseudo-spin quantization frame is given by the orthonormal unit vectors 
X, Y , Z , where Z is the quantization axis. In section 4.4 we set Z = A.z/d. 

• The i^^ ion's vibration around its equilibrium position is expressed in the 
local coordinate frame mf, see (12). For the Kitaev model we use the vectors 
given in (33): each vibrational direction (depending on which sublattice the ion is 
located) is indexed by, and associated with, one of the spin-space directions X, Y, 
Z, but this does not mean that the vibrational directions are parallel (or in any 
way related) to the spin-space axes. 

• The vectors connecting neighboring ions in the Kitaev honeycomb lattice Ax, 
Ay, A.Z are of length d and given in (32). They all lie in the plane of the lattice 
and do not form a 3D coordinate system. 
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